곡선 맞춤

물리학 실험이나 자연에 대한 관찰을 통해 얻어지는 데이터는 여러 다양한 이유로 오류가 포함되어 있어 정확하지 않은 경우가 많다. 예를 들어, 먼 은하까지의 거리 [eqn: d]와 그 은하가 우리로부터 멀어지는 속도 [eqn: v]를 측정했다고 해보자. 두 양이 서로 비례해서, 먼 은하일수록 더 빨리 멀어지고 있다는 것이 바로 천문학자 허블이 밝힌 유명한 허블의 법칙이다. 두 양 사이의 관계를 수식으로 적으면 [eqn: v`=`H`d]가 되는데, 실제의 관측 자료로 얻은 [eqn: d]와 [eqn: v]가 이 식을 따라 정확히 직선위에 놓이지는 못한다. 곡선맞춤(curve fitting) 방법을 이용하면, 관측 자료 [eqn: d]와 [eqn: v]가 [eqn: v`=`H`d]를 따른다는 가정 하에 허블 상수 [eqn: H]의 값으로 어떤 것이 가장 적절한지를 구할 수 있다.

일반적으로 두 양 사이에 [eqn: y`=`f(x)]의 관계가 있다고 하자. 실험이나 관측을 통해 얻어지는 [eqn: N]개의 자료는 [eqn: LEFT{(x _{1} ,y _{1} ),`(x _{2} ,y _{2} ),`(x _{3} ,y _{3} ), CDOTS ,`(x _{N-1} ,y _{N-1} ),`(x _{N} ,y _{N} ) right]의 형태로 적을 수 있다. 곡선맞춤은 주어진 데이터 자료를 이용해 가장 적절한 함수 꼴 [eqn: f(x)]를 정하는 것을 뜻한다. 실제 데이터와 가장 가까운 함수를 찾기 위해서는 실험값과 함숫값 사이의 차이를 가장 작게 하면 된다. 가장 널리 쓰이는 방법은 최소제곱법(least-square method)이라 불린다. ‘제곱’이라는 말이 들어있는 이유는 관측값과 함숫값, 둘 사이의 차이의 제곱을 모두 더한 양을 이용한다는 뜻이고, ‘최소’의 의미는 바로 이 양을 최소로 한다는 뜻이다. 실제 자료와 함숫값 사이의 차이의 제곱을 모두 더한 양을 오차(error)의 뜻으로 [eqn: E]라 적으면, [eqn: E = sum_{n=1}^N [y_n - f(x_n )]^2]의 형태로 적을 수 있다. 대부분의 상황에서는 함수 [eqn: f(x)]의 형태는 주어져있고, [eqn: E]를 최소로 하는 함수에 들어있는 조절변수의 값을 구하는 문제가 된다. 예를 들어 [eqn: f(x) = a+bx]와 같은 형태라면, 최소제곱법을 이용해 찾고자 하는 것은 [eqn: a]와 [eqn: b]의 값이다. 함수가 가진 이와 같은 조절변수가 [eqn: M]개라면 조절변수 전체를 하나의 벡터 [eqn: {vec{a}} =(a _{1} ,`a _{2} ,`a _{3} , CDOTS ,`a _{M} )]처럼 적을 수 있다. [eqn: E]의 최솟값은 모든 [eqn: a_i]에 대해서 [eqn: {partial E} over {partial a_i} = 0]를 만족한다. 가장 간단한 예로 [eqn: f(x)=bx]라면, [eqn: E= sum _{n=1} ^{N} (y _{n} -bx _{n} ) ^{2}]이므로 이를 미분하면,  [eqn: ` {Partial E} over {Partial b} =`-2 sum _{n=1} ^{N} (y _{n} -bx _{n} )x _{n} =0]이다. 즉, [eqn: b = {sum x_n y_n} over {sum x_n^2 }]를 얻게 된다. 위에서 소개한 허블 상수 [eqn: H]를 구하는 문제가 바로 이 경우다. 교재의 홈페이지에서 hubble.txt를 내려받아서 허블 상수를 구해보자. 파일 안의 각 줄의 첫 숫자는 Mpc의 단위로 적은 은하까지의 거리이며, 두 번째 숫자는 km/s의 단위로 적은 은하가 멀어지는 속도다. 참고로, 파일에 들어있는 데이터는 허블의 초기 논문의 값이라서 현재 우주론에서 추정하고 있는 값과는 큰 차이가 있다.

import numpy as np, matplotlib.pyplot as plt

data = np.loadtxt("hubble.txt")

x = data[:,0]

y = data[:,1]

Sxy = np.sum(x*y)

Sxx = np.sum(x*x)

H = Sxy/Sxx

print(H)

plt.plot(x,y,"x")

plt.plot(x,H*x)

plt.show()

위 프로그램에서는 식 [eqn: b = {sum x_n y_n} over {sum x_n^2 }]을 직접 적용해 허블 상수 [eqn: H]를 구했다. 허블의 원자료 데이터와 곡선 맞춤법을 적용해 얻어진 함수를 그래프로 함께 볼 수 있다. 허블 상수 [eqn: H]는 시간의 역수의 차원을 가진다. 현재 우리로부터 가장 먼 은하까지의 거리를 [eqn: R], 그 은하의 후퇴속도를 [eqn: V]라 하고, 그 은하가 우리와 같은 위치에 있었던 순간으로부터 경과한 시간을 [eqn: T]라 하자. 만약 [eqn: V]가 일정했다고 가정하면, [eqn: T=R/V]인데, 허블의 법칙 [eqn: v = Hd]에 의하면, 이 값이 다름 아닌 허블 상수의 역수가 된다. 위 프로그램에서 얻은 허블 상수를 가지고 우주의 나이를 추정하면 얼마가 될까 계산해 보라.

== 예제 1==

교재의 홈페이지에서 “fish.txt”파일을 내려 받으시오. 이 파일에는 우리나라 민물고기 한 종의 여러 개체의 길이(첫 컬럼, [eqn: L])와 무게(둘째 컬럼, [eqn: M])의 측정 자료가 들어있다. 무게와 길이는 [eqn: M=AL ^{b}]의 관계가 있고, [eqn: y = ln M], [eqn: x = ln L]로 바꿔 적으면 [eqn: y=a+bx]의 꼴이 된다. 오차 함수 [eqn: E]로부터 [eqn: Partial E/ Partial a=0,` Partial E/ Partial b=0]을 얻고 이를 연립해서 풀어 [eqn: E]가 최소가 되는 [eqn: a]와 [eqn: b] 각각의 표현식을 유도하시오. 내려 받은 파일에 들어있는 데이터를 이용해, [eqn: a]와 [eqn: b]의 값을 구하는 프로그램을 작성하시오. 자료 데이터와 얻어진 함수를 함께 그래프로 그리시오.

최소제곱법을 적용해 곡선 맞춤을 하는 위에 소개한 방법에서는 오차 함수 [eqn: E]가 최소가 될 때 그  미분값이 0이라는 것을 이용했다. 데이터를 기술하는 함수의 모양이 복잡하거나, 함수에 들어있는 조절 변수가 여럿일 때는, [eqn: partial E /partial a_i = 0]을 모든 [eqn: a_i (i=1, 2, cdots, M)]에 대해서 연립해서 풀기가 쉽지 않은 경우가 많다. 이럴 때는, [eqn: vec a]의 현재 값에서 출발해서 [eqn: E]를 조금씩 한 단계씩 줄여나가는 방법을 많이 쓴다. 미적분학에서 배운 것을 돌이켜 보자. 스칼라 마당(scalar field) [eqn: phi(vec x )]이 주어져 있을 때 [eqn: vec x + d vec x]에서의 값과 [eqn: vec x]에서의 값의 차이는 [eqn: d phi = phi(vec x + d vec x ) - phi(vec x )= del phi cdot d vec x]로  경사(gradient, 기울기라고도 부름) [eqn: Del phi]를 이용해 적을 수 있다. 이 식에 의혐, [eqn: d vec x]의 방향을 [eqn: del phi]와 같은 방향으로 고르면, [eqn: phi(vec x )]가 가장 빨리 증가하는 방향이 된다. 오차 함수 [eqn: E]를 줄이는 수치 계산 방법인 경사 하강법(gradient descent)이 바로 이 식을 이용한다. 물론, [eqn: E]를 늘이는 것이 아니라 줄이는 것이 목표이므로, [eqn: - del E]의 방향을 택하면 된다. 자, 조절 변수들([eqn: a_1 , a_2 , cdots , a_M])이 만드는 벡터공간에서 현재의 위치가 [eqn: {vec{a}}]라면, 가장 빨리 [eqn: E]의 값을 줄이는 방향은 바로 [eqn: -Del E]의 방향이다. 즉, 적당히 작은 상수 [eqn: c]를 택하고 [eqn: d {vec{a}} =`- c NABLA E]로 하면, 각각의 조절 변수 [eqn: a _{i}]에 대해 [eqn: a _{i} `` -> ```a _{i} -`c` {Partial E} over {Partial a _{i}}]의 과정을 반복하해 오차 함수 [eqn: E]가 최소인 점을 찾아갈 수 있게 된다. 오차함수 [eqn: E = sum_{n=1}^N [y_n - f(x_n )]^2]을 조절변수의 함수로 생각해 적으면 [eqn: E( {vec{a}} )= sum _{n=1} ^{N} [y _{n} -f(x _{n} `;` {vec{a}} )] ^{2}]의 꼴이므로 [eqn: {partial E} over{ partial a_i }= - 2 sum [ y_n - f(x_n ; vec a )] {partial f(x_n ; vec a )} over {partial a_i }]이다. 자, 경사하강법을 따라 위에서 살펴본 허블 상수의 문제를 구현해 보자. [eqn: f(x;H)=Hx]로 적을 수 있으므로, [eqn: {Partial E} over {Partial H} =-2 sum _{} ^{} (y _{n} - H x _{n} ) x_n]이다. 즉, 현재의 [eqn: H]값을 다음 단계에서는 [eqn: H + c sum _{} ^{} (y _{n} -Hx _{n} )x _{n}]로 바꾸는 과정을 반복하면 된다. 이때 [eqn: c]는 충분히 작은 값을 택해서 각 단계에서 [eqn: H]가 변하는 정도가 너무 크지 않도록 하는 것이 좋다.

import numpy as np, matplotlib.pyplot as plt

eps = 1.0E-5; c = 0.01; H = 10.0

data = np.loadtxt("hubble.txt")

x = data[:,0]; y = data[:,1]

Sxy = np.sum(x*y); Sxx = np.sum(x*x)

dH = c*(Sxy - H*Sxx)

while ( np.abs(dH) > eps ):

H += dH

dH = c*(Sxy - H*Sxx)

print(H)

plt.plot(x,y,"x")

plt.plot(x,H*x)

== 예제 2==

예제 1의 [eqn: a, b]를 경사 하강법을 이용해 계산하는 프로그램을 작성하시오. 자료 데이터와 얻어진 함수를 함께 그래프로 그리시오.

== 과제 1==

4장의 과제 1에서 “mers.txt”파일의 데이터와 미분방정식은 [eqn: dN over dt = r N left( 1 - N over K right)]의 결과를 비교한 바 있다.  [eqn: N(t=0) = 1]을 처음 조건으로 하고 이 미분방정식의 해석적인 해 [eqn: N(t;`r,K)]를 구하시오. (프로그램과는 별도로 풀이를 적어 제출하시오). [eqn: r]과 [eqn: K]의 값을 경사 하강법을 이용해 구하는 프로그램을 작성하시오. 데이터와 최종 곡선 맞춤 결과를 그래프로 함께 표시하시오.

위에서 경사 하강법을 이용해 곡선 맞춤을 하는 방법을 소개했다. 이를 위해서는 [eqn: {partial E} over{ partial a_i }= - 2 sum [ y_n - f(x_n ; vec a )] {partial f(x_n ; vec a )} over {partial a_i }]를 계산해야 했다. [eqn: f(x; vec a )]의 편미분 [eqn: {Partial f(x _{n} ; {vec{a}} )} over {Partial a _{i}}]을 해석적으로 구하고, 이를 프로그램에서 직접 이용했다. 다른 방법도 있다. 프로그램 안에서 수치적인 방법으로 미분값의 어림값을 구해도 된다. 방법은 간단하다. 현재의 [eqn: vec a]에서 [eqn: a_i]를 [eqn: a_i + da_i]로 변화시킨 후, [eqn: {Partial E} over {Partial a _{i}} APPROX {E(a _{1} ,`a _{2} , CDOTS ,`a _{i} +da _{i} ,`a _{i+1} , CDOTS ,`a _{M} )- E(a _{1} ,`a _{2} , CDOTS ,`a _{i} ,`a _{i+1} , CDOTS ,`a _{M} )} overda_i]을 이용하면 된다. 허블 상수를 구하는 위의 프로그램을 수치 미분을 적용해 다시 적어보자.

import numpy as np, matplotlib.pyplot as plt

def f(x,H):

return(H*x)

def E(x, y, H):

return np.sum( (y - f(x,H))**2 )

eps = 1.0E-5

c = 0.01; dH = 0.001

data = np.loadtxt("hubble.txt")

x = data[:,0]

y = data[:,1]

H = 10.0; DeltaH=1

while ( np.abs(DeltaH) > eps ):

dEdH = (E(x,y,H+dH) - E(x,y,H))/dH

DeltaH = -c*dEdH

H += DeltaH

print(H)

plt.plot(x,y,"x")

plt.plot(x,H*x)

plt.show()

위 프로그램에서 수치 미분을 위해 충분히 작은 H의 변화량 dH를 정의해 이용했다. 또한, 한 단계에서의 H의 변화량 DeltaH가 eps로 정의된 값보다 작아지면 계산을 멈추도록 했다. DeltaH는 while loop안에서 계산되므로 처음 값으로 적당히 큰 값 1을 부여했다.

== 예제 3==

예제 2의 프로그램을 수치 미분을 이용하도록 바꿔 작성하시오.

== 과제 2==

과제 1의 프로그램을 수치 미분을 이용하도록 바꿔 작성하시오.